%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code to recover the hzd rates in 3-state model of the labor market
% R. Barnichon and G Mesters, 03/07/2017
% 5 steps:
% 1) input measured transition probas
% 2) correct for 1994 redesign
% 3) correct for time aggregation bias
% 4) correct for magin error
% 5) put all hzd rates together in one mat file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear all

addpath ('var_bvar')

% names of different demographic groups
a=['w1624'; 'w2534'; 'w3544'; 'w4554'; 'w5564'; 'm1624'; 'm2534'; 'm3544'; 'm4554'; 'm5564'; '_6585'];


lambda=[];

for index=1:length(a)
    
    %load raw CPS transition data
    data=xlsread(['Flows\Flows_',num2str(a(index,:)),'.xlsx']);
    data=data(1:end-1,:);
    
    neu=data(:,1);
    nei=data(:,2);
    niu=data(:,3);
    nie=data(:,4);
    nue=data(:,5);
    nui=data(:,6);
    
    % order of shares:
    n=[neu nei niu nie nue nui];
    
    %save origina\l series:
    n_init=n;
    
    E=data(:,7);
    U=data(:,8);
    I=data(:,9);
    
    % je vire les elements vides et j'interpole
    for j=1:size(n,2)
        for i=2:length(n)-1
            if n(i,j)==0 | isnan(n(i,j))==1
                if (n(i-1,j)~=0 && n(i+1,j)~=0) & (isnan(n(i-1,j))~=1 && isnan(n(i+1,j))~=1)
                    n(i,j)=(n(i-1,j)+n(i+1,j))/2;
                else
                    n(i,j)=n(i-2,j);
                end
            end
        end
    end
    
    neu=n(:,1);
    nei=n(:,2);
    niu=n(:,3);
    nie=n(:,4);
    nue=n(:,5);
    nui=n(:,6);
    
    
    % XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    % CPS REDESIGN (94) CORRECTION
    % XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    
    %1994m1
    date=(1994-1976)*12+1-1; %because we start in 76m2
    
    Y=1994;M=1;
    actualend=(Y-1976)*12 +M-1; %ie, first point of forecast
    
    % XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    % Use VAR to forecast
    % XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    name=cd;
    addpath ([name,'../var_bvar'])
    Je=[];
    Je=[log(n)];
    Tf=42;
    
    % We only observe data until actualend-1 for the flows
    Je=Je(1:actualend-1,:);
    
    %  Only use the last 10 years of data to estimate the VAR ("rolling estimation window")
    Jet=Je(actualend-121:actualend-1,:);
    
    % NUMBER OF LAGS: I USE 6
    Jf=varf(Jet,6,Tf,length(Jet)+1);
    Jfor=[Je(1:actualend-1,:)' Jf']';
    
    % figure, plot(Jfor);
    % hold on, plot([log(jeu) log(jue) log(jei) log(jie) log(jiu) log(jui) ])
    % hold on, plot(Je(1:end-Tf,:));
    
    lNfor=(Jfor);
    
    forlength=12;
    
    lNfor_avg94=mean(lNfor(actualend:actualend+forlength,:));
    ln_avg94=mean(log(n(actualend:actualend+forlength,:)));
    
    adj_level=lNfor_avg94-ln_avg94;
    
    ln_alt=log(n);
    
    % Adjust data pre-1994:
    ln_alt(1:actualend,:)=log(n(1:actualend,:))-kron(adj_level,ones(length(n(1:actualend,:)),1));
    
    n_alt=exp(ln_alt);
    %     figure, plot([n n_alt])
    
    n=n_alt;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TAB Adjustment
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    for i=2:length(n)
        neu=n(i,1);
        nei=n(i,2);
        niu=n(i,3);
        nie=n(i,4);
        nue=n(i,5);
        nui=n(i,6);
        
        % use method of EHS
        % Write transition matrix: (ordering E, U, I)
        P=[1-neu-nei nue nie;...
            neu 1-nue-nui niu; ...
            nei nui 1-nie-niu;];
        
        [V D]=eig(P);
        lnD=log(diag(D));
        D=diag(lnD);
        F=real(V*D*inv(V));
        
        leu(i,index)=F(2,1);
        lei(i,index)=F(3,1);
        
        lue(i,index)=F(1,2);
        lie(i,index)=F(1,3);
        
        lui(i,index)=F(3,2);
        liu(i,index)=F(2,3);
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Take Q averages of M data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %I take qtly averages but I ignore the first three months of 1976 (feb and March) to start in 1976Q2
    leu_q=[];
    leu_q=mq(leu(4:end,index));
    leu_q=[(leu(2,index)+leu(3,index))/2; leu_q];
    
    lei_q=mq(lei(4:end,index));
    lei_q=[(lei(2,index)+lei(3,index))/2; lei_q];
    
    liu_q=mq(liu(4:end,index));
    liu_q=[(liu(2,index)+liu(3,index))/2; liu_q];
    
    lie_q=mq(lie(4:end,index));
    lie_q=[(lie(2,index)+lie(3,index))/2; lie_q];
    
    lue_q=mq(lue(4:end,index));
    lue_q=[(lue(2,index)+lue(3,index))/2; lue_q];
    
    lui_q=mq(lui(4:end,index));
    lui_q=[(lui(2,index)+lui(3,index))/2; lui_q];
    
    
    l= [leu(:,index) lei(:,index) liu(:,index) lie(:,index) lue(:,index) lui(:,index)];
    n=1-exp(-l);
    
    
    % XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    % MARGIN ERROR CORRECTION
    % XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    
    n=1-exp(-l);
    
    %Set 1976m1 to NaN
    n(1,:)=nan;
    n_alt(1,:)=nan;
    
    %Post 1976m1:
    for i=2:length(n)
        neu=n(i,1);
        nei=n(i,2);
        niu=n(i,3);
        nie=n(i,4);
        nue=n(i,5);
        nui=n(i,6);
        El=E(i-1);
        Ul=U(i-1);
        Il=I(i-1);
        
        % I use the weighting matrix of Elsby et al. (EHS, 2013)
        W=inv([neu*(1-neu)/E(i-1) -neu*nei/E(i-1) 0 0 0 0; ...
            -neu*nei/E(i-1) nei*(1-nei)/E(i-1) 0 0 0 0; ...
            0 0 nue*(1-nue)/U(i-1) -nue*nui/U(i-1) 0 0; ...
            0 0 -nue*nui/U(i-1) nui*(1-nui)/U(i-1)  0 0; ...
            0 0 0 0 nie*(1-nie)/I(i-1) -nie*niu/I(i-1); ...
            0 0 0 0 -nie*niu/I(i-1) niu*(1-niu)/I(i-1) ]);
        
        X_1=[-E(i-1) -E(i-1) U(i-1) 0 I(i-1) 0; ...
            E(i-1) 0 -U(i-1) -U(i-1) 0 I(i-1)];
        
        p_hat=[neu;nei;nue;nui;nie;niu];
        
        ds=[E(i)-E(i-1);U(i)-U(i-1);];
        
        %Lagrange multiplier:
        mu=inv(X_1*inv(W)*X_1')*(ds-X_1*p_hat);
        % corrected transition proba:
        p=p_hat+inv(W)*X_1'*mu;
        
        % Check that flows are consistent with the stocks:
        X_1*p-ds
        if abs(X_1*p-ds)>.000000001
            i
            pause
        end
        n_alt(i,:)=[p(1) p(2) p(6) p(5) p(3) p(4)];
        
        %Get rid of probas above one
        for j=1:size(n_alt,2)
            if n_alt(i,j)>1
                n_alt(i,j)=n_alt(i-1,j);
            end
        end
        
    end
    
    n=n_alt;
    
    l=-log(1-n);
    
    leu(:,index)=n(:,1);
    lei(:,index)=n(:,2);
    liu(:,index)=n(:,3);
    lie(:,index)=n(:,4);
    lue(:,index)=n(:,5);
    lui(:,index)=n(:,6);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %Get monthly E, U and I:
    Em=(E);
    Um=(U);
    Im=(I);
    
    %Get quarterly E, U and I:
    Eq=mq(E);
    Uq=mq(U);
    Iq=mq(I);
    
    %transition proba:
    p=1-exp(-l);
    peu=p(:,1);
    pei=p(:,2);
    piu=p(:,3);
    pie=p(:,4);
    pue=p(:,5);
    pui=p(:,6);
        
    size(E)
    size(l)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Take Q averages of M data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    leu_q=[];
    leu_q=mq(leu(4:end,index));
    leu_q=[(leu(2,index)+leu(3,index))/2; leu_q];
    
    lei_q=mq(lei(4:end,index));
    lei_q=[(lei(2,index)+lei(3,index))/2; lei_q];
    
    liu_q=mq(liu(4:end,index));
    liu_q=[(liu(2,index)+liu(3,index))/2; liu_q];
    
    lie_q=mq(lie(4:end,index));
    lie_q=[(lie(2,index)+lie(3,index))/2; lie_q];
    
    lue_q=mq(lue(4:end,index));
    lue_q=[(lue(2,index)+lue(3,index))/2; lue_q];
    
    lui_q=mq(lui(4:end,index));
    lui_q=[(lui(2,index)+lui(3,index))/2; lui_q];
    
    l_q= [leu_q lei_q liu_q lie_q lue_q lui_q];
    figure, plot(l_q)
    name=['Flows\J_2016',num2str(a(index,:)),'.mat'];
    save (name, 'l','l_q','Em','Um','Im')
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Put all monthly hzd rates in one mat file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

JEU=[];
JEI=[];
JIU=[];
JIE=[];
JUE=[];
JUI=[];
E=[];
U=[];
I=[];

a=['w1624'; 'w2534'; 'w3544'; 'w4554'; 'w5564'; 'm1624'; 'm2534'; 'm3544'; 'm4554'; 'm5564'; '_6585'];

for index=1:length(a)
    name=['Flows\J_2016',num2str(a(index,:)),'.mat'];
    load (name, 'l','Em','Um','Im');
    
    JEU=[JEU l(:,1)];
    JEI=[JEI l(:,2)];
    JIU=[JIU l(:,3)];
    JIE=[JIE l(:,4)];
    JUE=[JUE l(:,5)];
    JUI=[JUI l(:,6)];
    E=[E Em];
    U=[U Um];
    I=[I Im];
    
end

% Get rid of negative values
for j=1:size(JEI,2)
    for i=2:size(JEI,1)
        if JEI(i,j)<=0, JEI(i,j)=JEI(i-1,j); end
        if JEU(i,j)<=0, JEU(i,j)=JEU(i-1,j); end
        if JIE(i,j)<=0, JIE(i,j)=JIE(i-1,j); end
        if JIU(i,j)<=0, JIU(i,j)=JIU(i-1,j); end
        if JUI(i,j)<=0, JUI(i,j)=JUI(i-1,j); end
        if JUE(i,j)<=0, JUE(i,j)=JUE(i-1,j); end
    end
end


save('Hzd.mat','JEU','JEI','JUE','JUI','JIE','JIU','E','U','I')

